Noise enhanced persistence in a biochemical regulatory network with feedback control 
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We find that discrete noise of inhibiting (signal) molecules can greatly delay the extinction of 
plasmids in a plasmid replication system: a prototypical biochemical regulatory network. We calcu- 
late the probability distribution of the metastable state of the plasmids and show on this example 
that the reaction rate equations may fail in predicting the average number of regulated molecules 
even when this number is large, and the time is much shorter than the mean extinction time. 
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Many molecular species that control genetic regula- 
tory networks are present in low concentrations. The 
resulting fluctuations in reaction rates may cause large 
random variations in the instantaneous intracellular con- 
centrations of molecular species which, in their turn, 
may have important consequences in biological function- 
ing. This and related topics have attracted much re- 
cent interest from the biology and physics communities 
P, 0, 0, 0, 0, 0, 0, 0] • Intracellular processes are often reg- 
ulated via negative feedback by signal molecules. It was 
assumed in the past that noise in the signal component 
would randomize control of the regulated component. 
More recently, it has been shown that this noise may 
actually enhance the robustness of the regulated compo- 
nent, bringing the variation of its probability distribution 
below the Poissonian limit [J]. Here we report a previ- 
ously unexplored dramatic impact the noise can have on 
the persistence of the regulated component in systems 
with negative feedback control. Following Paulsson et 
al. we will consider a minimal two-component copy 
number control (CNC) model that, on the one hand, in- 
cludes standard intracellular processes and, on the other 
hand, provides an adequate description to CNC of bac- 
terial plasmids. Plasmids are extra-chromosomal DNA 
molecules (typically, circular and double-stranded) that 
are capable of autonomous replication. They undergo in- 
tracellular dynamics of the birth-death type with decay 
(mostly dilution by cell division) and autocatalytic pro- 
duction inhibited by signal molecules. If there is no pen- 
etration of new plasmids into the cell, a rare sequence 
of multiple decay events will ultimately drive the plas- 
mid population to extinction. It may even cause the 
death of the cell if the plasmid contains a vital gene. 
We will combine analytical and numerical approaches to 
show that noise in the number of signal molecules can 
greatly delay the plasmid extinction. We will also cal- 
culate analytically the probability distribution function 
(PDF) of the metastable state of the regulated molecules 
and show that widely used deterministic reaction rate 
equations (RRE) may fail in predicting the average num- 
ber of plasmids even when this number is large. These 
remarkable effects do not require an unusual molecular 
distribution and occur due to rare events when the num- 



ber of signal molecules is very small. 

Model. Consider a double negative-positive feedback 
loop with plasmids denoted by X and signal molecules 
denoted by S. The plasmids promote the production of 
the signal molecules, whereas the signal molecules inhibit 
the autocatalytic production of the plasmids. The RRE 
for the average concentrations of the two species are [4j: 



X = X^(S/A)-X, 
S = aX-/3S, 



(1) 



where ^(S/A) is a nonlinear and monotone decreasing 
function of S, ^(O) > 1, the parameter A specifies the 
inhibition strength of the signal molecules S, and time 
and the rates are rescaled by the decay rate of the plas- 
mids. Equations (p} have an attracting fixed point (X, S) 
[where X = S/rj, ty(S/A) = 1, and r] — a/ [3] and an un- 
stable fixed point (X ,Sq) = (0,0). According to the 
RRE, the system would stay in the (X, S) state for- 
ever. The underlying stochastic process, however, be- 
haves quite differently. A large enough fluctuation ulti- 
mately depletes the plasmid population. The state with 
no plasmids is an absorbing state, as the probability of 
escape from it is zero. Therefore, the (Xq,Sq) state is 
actually stable, whereas the stable fixed point (X, S) of 
the deterministic model is metastable. The mean extinc- 
tion time (MET): the mean time it takes this stochastic 
process to reach the absorbing state is expected to be ex- 
ponentially long in the (presumably large) average num- 
ber of plasmids in the metastable state, see e.g., Refs. 

To account for the stochastic effects, consider a chem- 
ical master equation (CME) that describes the evolution 
of the probability P m , n {t) of having, at time t, m plas- 
mids and n S-molecules. For m, n > 1 the CME is 



+ am{E- 1 - l)P m<n + - l)nP m 



(2) 



where = f(n + j) and g m>n = m^(n/A) [lj. 

Let us denote by P n \ m the probability of having, at time 
t, n S'-molecules conditioned on having m plasmids, and 
by 7T m the probability of having m plasmids regardless of 
the number of S'-molecules. We substitute the identity 
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FIG. 1: The phase plane X, S of the reaction rate equations 
Q for V(S/A) = kexp{-S/A), = 150, a = 200, A = 4, and 
fc = 14. The circle denotes the fixed point (X, S). 



Pm,n = P^m^m into Eq. © and sum over all n. The 



result is 



7T m = - l)g m 7r m + - l)rmr„ 



(3) 



where g m = J2n=o P n\m(t)9m,n is the production rate 
of the plasmids averaged over the (yet unknown) condi- 
tional distribution P n \ m (t). 

Further analytical progress is only possible in some lim- 
its. Following Paulsson et al. [3], we assume that the 
S'-dynamics is much faster than the X-dynamics. At the 
level of the RRE, S adjusts rapidly to the current value 
of X. Then S{t) ~ T}X{t) holds, while X(t) and S(t) flow 
relatively slowly towards the fixed point (X, S) according 
to the reduced equation X ~ X^(nX/A) - X. We will 
perform further calculations in two particular examples: 
exponential and hyperbolic inhibition. 

Exponential inhibition. Here ^(S/A) — k exp(— S/A), 
k > 1 (we assume that k is not too close to 1), and 
(X,S) = (Alnk/r), Alnfc). The time scale of the fast 
dynamics is ~ 1//3, the time scale of the slow dynamics 
is ~ 1/ In k [see Fig.Q] , so the time scale separation occurs 
when In k -C /3. 

At the level of the CME we can perform adiabatic elim- 
ination of the fast dynamics in the variable P n \m{t) by 
assuming that the S-population rapidly adjusts to a Pois- 
son distribution PjQ = e _T,m '(r)m) n / 'n\ about the cur- 
rent value of the mean, r)m(t). Now the effective stochas- 
tic rate g m in Eq. ([3|) can be easily calculated 0|: 



yp^c 



kme- rT > m , r = 1 - e 



-1/A 



(4) 



n=0 



This procedure reduces the two-species problem to an 
effective one-species problem: a single-step birth-death 
process with the birth rate g m and death rate \i m = m. 
In the most interesting case of (m) ^> 1, there are two 
widely different time scales in this process. The first, 
short time scale is the relaxation time to the metastablc 
state. The second, exponentially long, time scale is the 
life time of the metastable state, or the extinction time, 
see below. At intermediate times one observes a qua- 
sistationary distribution (QSD) q m (t) — Tm/[1 — TR) Wi- 



the PDF of having, at time t, m plasmids conditioned 
on their non-extinction, see e.g. |13j |. When gi/^i = 
ke~ rv 3> f, that is lnfc ^> rr], the probability flux to the 
zero state m = is negligible, and q m (t) can be approx- 
imated by putting, in Eq. ([3]), ~k m (t) = for all m and 
assuming [i\ = 11411. In this way one obtains a recursion 
relation for q m [lfj, [l5| which yields the QSD: 

„ — (1/2) r'qm(m—l) h.m—1 

—, (5) 
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while gi can be found from the normalization 
Sm=i 1™ ~ Assuming rr\ <C _L_we can replace the 
normalization sum by an integral [1 61 ] and, by the saddle 
point method, obtain 



cil 1 = r 



e 2 



\fk In k 

Now, <?f 1 is nothing but the MET r [13], and Eq. © 
yields an accurate approximation for it. The same result 
follows from an exact expression for the MET [13] • 

Let us calculate for comparison the MET for the "semi- 
deterministic" (SD) case: when S — i]X is a prescribed 
deterministic quantity. The SD rate = kme~ vm ^ A is 
obtained by putting n = r\m. A similar calculation, for 
r\j A <C f , yields 



(7) 



V kA In k 

How do the fully stochastic ^ and SD ([7]) results for the 
METs compare? Consider their ratio 

"(lnfc) 2 



R : 



'rA exp 



2rrj 



-(1-rA) 



(8) 



The strongest effect is observed for (lnfc) 2 ^S> r\. In 
this case, and for A n, we obtain R ^S> 1: the dis- 
crete noise of the S'-molecules greatly (exponentially) 
delays the plasmid extinction. Note that in this pa- 
rameter regime R is a monotone decreasing function of 
A. However, even for A — > oo the effect is strong, as 

R _ e (Infc) a /(^) » 1. 

Using Eq. ^ for r, we can determine the extinc- 
tion probability 7To(t): the probability that extinction 
occurs until time t, see e.g. Ref. [11|. Also, by con- 
servation of probability, we can restore the exponentially 
slow time-dependence of the PDF of the metastable state, 
7Tm>o(£) — q m exp(— t/r) [where r is given by Eq. ©]. 
We obtain 



T0(*) 
7Tm>o(i) 



l- e -*/ r , 
fc m-l/2 lnfc 



(9) 



Using the PDF ([9]) , we can calculate the (slowly decaying 
in time) average number of the plasmids: 

E /lnfc 1" 
rmr m (t) ~ h 

m=0 v ' 



(m(i)) 



-t/r 



(10) 



FIG. 2: (Color online) The PDF ir m (t) found by solving nu- 
merically a truncated CME ((2| for the exponential inhibition 
with k = 13, A = 4 and a — /3 = 400. The dashed line shows 
the initial distribution: a Kroenecker delta at m = n = 18. 
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FIG. 3: (Color online) (a) The PDF J9} of the metastable 
state for the exponential inhibition (solid line) and numerical 
solution of the CME ([2]) (circles) for f<T. The parameters 
are k = 13, A = 3 and a = (3 = 500. (b) The MET versus 
A~ l , the rest of parameters the same as in (a). Solid line: Eq. 
((6]), circles: numerical solutions for the fully stochastic case, 
dashed line: Eq. ([7]), squares: numerical solution for the SD 
case, see text for details. The ratio of the fully stochastic and 
SD METs increases with the inhibition strength 1/A. 



where we have again assumed rr\ <C 1. For A < 1, (m(t)} 
strongly deviates from the RRE prediction X = A \nk/n, 
even at t <C r, as was previously observed numerically [4j. 
Note that non-gaussianity of the PDF © appears only 
in the pre-exponent. 

To test our analytical results, we solved numerically a 
truncated CME @. The numerically found PDF of the 
plasmids 7r ro (i) exhibits a slow decay of the metastable 
state and a simultaneous growth of the extinction proba- 
bility in time, see Fig. [2] Figure [3] compares our analyti- 
cal and numerical results for 7r m (f ) . In addition, we com- 
pare there the analytical result ^ for the MET with the 
numerical result r num = -*/ln[l - E,T=o P o,»(i)] (that 
approaches a constant after a transient), and also T sd 
from Eq. ([7]) with the result of a numerical solution of 
Eq. §5§ with the SD rate Qm- Very good agreement is 
observed for all quantities [18j |. 

Hyperbolic inhibition. Our second example employs 
the widely used hyperbolic, or Michaelis-Menten, inhibi- 
tion model [ll]. Here V(S/A) = k/(l + S/A), k > 1 (and 



FIG. 4: (Color online) (a) The PDF of the metastable 
state for the hyperbolic inhibition (solid line) and numerical 
solution of the CME ([2]) (circles) for t<r. The parameters 
are k = 15, A = 1, a = 350, and /3 = 700. (b) The MET ver- 
sus A -1 , the rest of parameters the same as in (a). Solid line: 
Eq. (|12|1 . circles: numerical solutions for the fully stochastic 
case, dashed line: Eq. (|12|l with the SD g m , squares: numeri- 
cal solution for the SD case, see text for details. The ratio of 
the fully stochastic and SD METs grows with the inhibition 
strength 1/A. (c) Equation for the MET (solid line) and 
numerical solutions of the CME @ (circles) at different n for 
k = 10, A = 10~ 3 and f3 = 2 x 10 3 . 



not too close to 1), and (X, S) = [(k - l)A/r), (k - I) A]. 
The time scale separation occurs at j3 » 1. At the level 
of the CME (J2) we again assume a rapid adjustment of 
the 5-species to the m-dependent Poisson distribution. 
The effective stochastic rate g m in Eq. ([3]) is 



yp (P) gn 

n=0 



kme' 11 " 1 1 F 1 (A,A+ 1,77m), (11) 



where iF±(a,b, z) is the Kummer confluent hypergeo- 
metric function [2(|. Using this effective rate, Pauls- 
son and Ehrenberg Q calculated the QSD numerically. 
We have found it analytically from the recursion rela- 
tion [l(| fl5j |. by assuming gi » /Ji = 1. The result is 
Im/li — (1/m!) n5=i 9j- Again, q^ 1 = r can be found 
by normalizing the QSD to unity. Therefore, the PDF of 
having m plasmids at time t, and the MET, are 



-t/r 



m — 1 



00 m—1 



lift- ^e^lu- (12) 

i=i 



mi 

m=l j=l 



Comparisons between these predictions for the fully 
stochastic and SD cases [with truncated sums in Eq. (fl~2"j) ] 
and numerical solutions of the truncated CMEs (J2]) and 
([3]), respectively, are shown in Fig. 0^ and b, and very 
good agreement is observed [Tij . 

The extreme case of very strong inhibition, A <C 
ln/e/fc, can be further simplified. It can be checked a 
posteriori that here, for all m that contribute to the 
normalization of 7r m , and hence to the MET, the ef- 
fective rate (fTTj) is well approximated by the first term: 
g m ~ kme^ r,m . This rate formally coincides with that 
given by Eq. Q for the exponential inhibition, where 
one must put r = 1. Therefore, the most interesting 
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case of strong inhibition A <C Ink/k (when the stabi- 
lizing effect of noise in the S-molecules on the plasmid 
fluctuations and persistence is the largest) is also the 
simplest. Furthermore, the exponential inhibition model 
formally describes the strong hyperbolic inhibition limit. 
By additionally assuming that A <C r)/k <C lnfc/fc, one 
can show after some algebra that the QSD and n m (t) 
from Eq. (fT2|) reduce to Eqs. and ([9|), respectively 
(with r = 1). The slowly-decaying average of this PDF 
[Eq. (fT0|) with r = 1] again strongly deviates, already at 
t < r, from the RRE prediction X = (k — l)A/r). The 
corresponding MET [Eq. ([6]) with r = 1] can again be 
compared with the SD MET, obtained by using the SD 
rate g^ = fcm(l + rjm/A)^ 1 which assumes n — T]m. 
For A <C rj/k one has gf d <C 1, so T s d = Oil), as the 
decay dominates over the replication. In contrast, the 
asymptotic result for the stochastic MET, for rj <C 1, 



(13) 



Vic In A: 

is exponentially large. Therefore, the noise in the number 
of the S'-molecules again causes, at A <C n/k, exponential 
enhancement of the persistence of the plasmids. Equa- 
tion (|13p compares well with numerics, see Fig. |4j;. 

Discussion. We have shown, in a simple CNC model, 
that intrinsic discrete noise of the signal molecules 
can greatly increase the average number of regulated 
molecules and therefore enhance the persistence of the 
regulated component. Although we assumed that P n \ m 
is Poisson distributed, we expect these findings to hold, 
for sufficiently strong inhibition, for other signal molecule 
kinetics as well. What is the mechanism behind the noise 
enhanced persistence? The autocatalytic production rate 
of the plasmids is largest at S — 0, therefore rare events 
of having a very small number of S'-molecules strongly 
dominate the effective stochastic growth rate g m , see e.g. 
Eq. (HJ). As a result, the average number of plasmids in 
the metastable state greatly increases, and this enhances 
the plasmid persistence. As the mode and the average for 
the plasmid PDF 7r m coincide, this mechanism of failure 
of the RRE is different from that discussed previously 0] • 

The noise-enhanced persistence, that we predict here, 
should be observable in experiment, in vitro and in vivo, 
due to recent advances in single-molecule signal mea- 
surements. Finally, the effect is not system-specific, and 
should appear in a host of other birth-death-type systems 
where negative feedback is at work. 

We thank Ari Meerson for a useful discussion. Our 
work was supported by the Israel Science Foundation. 
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